Average Bubble Profile¶
In [1]:
import numpy as np
import math
import matplotlib
#matplotlib.rcParams.update({'font.size': 12})
from matplotlib import gridspec
import matplotlib.pyplot as plt
from collections import deque
import scipy as scp
import scipy.optimize as sco
import scipy.signal as scs
import scipy.special as ssp
import scipy.integrate as sci
import scipy.interpolate as intp
import statistics as stat
import random
from functools import partial
from operator import eq
from itertools import zip_longest, compress, count, islice, groupby, cycle
#from labellines import labelLine, labelLines
import os
import scipy.ndimage
from scipy.ndimage import gaussian_filter, gaussian_filter1d
In [2]:
np_load_old = np.load
np.load = lambda *a,**k: np_load_old(*a, allow_pickle=True, **k)
Simulation Parameters¶
In [3]:
nLat = 1024
lSim = 0
nSims = 300
filter_size = 0.
temp = 0.
phi0 = 2.5
lamb = 1.4
nCols = 2
nu = 2.*10**(-3)
V0 = 4.*nu
m2eff = lambda lam: V0 * (- 1. + lam**2)
lenLat = 250. / np.sqrt(V0); print('lenLat = ', lenLat)
phi_initial = np.pi
mask = 4*phi_initial
#normal = [phi_initial, np.mean([sim[1][0] for sim in all_data]), np.mean([sim[2][0] for sim in all_data]), V(phi_initial), np.mean([0.5*sim[1][0]**2 for sim in all_data])]
#normal = [phi_initial, 0, 0, V(phi_initial), 0]
normal = [phi_initial, 0]
nyq = nLat//2+1; spec = nyq; dx = lenLat/nLat; dk = 2.*np.pi/lenLat; dtout=dx; alpha=8.; print('dx, dk, spec ', dx, dk, spec)
light_cone = int(dtout/dx); print('light_cone = ', light_cone)
V = lambda phi, lamb: V0 * (-np.cos(phi)+0.5*lamb**2.*np.sin(phi)**2.)
dV = lambda phi, lamb: V0 * (np.sin(phi)+0.5*lamb**2.*np.sin(2.*phi))
far_right_phi_at_V_max = sco.minimize_scalar(lambda x: -V(x, lamb), bounds=[2*np.pi, 3*np.pi], method='bounded')
right_phi_at_V_max = sco.minimize_scalar(lambda x: -V(x, lamb), bounds=[np.pi, 2*np.pi], method='bounded')
left_phi_at_V_max = sco.minimize_scalar(lambda x: -V(x, lamb), bounds=[0, np.pi], method='bounded')
phi_upper_bound = sco.fsolve(lambda x: V(x, lamb) - V(phi_initial, lamb), 5)[0]; print(phi_upper_bound)
phi_lower_bound = sco.fsolve(lambda x: V(x, lamb) - V(phi_initial, lamb), 1)[0]; print(phi_lower_bound)
plt.plot(np.linspace(0, 2*np.pi, 100), V(np.linspace(0, 2*np.pi, 100), lamb))
plt.plot(phi_upper_bound, V(phi_upper_bound, lamb), 'ro', phi_lower_bound, V(phi_lower_bound, lamb), 'go')
plt.plot(right_phi_at_V_max.x, V(right_phi_at_V_max.x, lamb), 'bo', left_phi_at_V_max.x, V(left_phi_at_V_max.x, lamb), 'ko')
Out[3]:
In [4]:
titles = [r'$\phi(x)$', r'$\partial_t \phi(x)$', r'$|\nabla \phi(x)|^2$', r'$V(\phi(x))$']
plots_file = '/home/dpirvu/project/thermal_bubbles/plots/'
pickle_file = '/home/dpirvu/project/pickle_location/thermal_bubbles/'
sim_location = lambda nL, tem, phi, lam, sim: '/gpfs/dpirvu/thermal_bubbles/bubbles_x'+str(nL)+'_temp{:.4f}'.format(tem)+'_phi0{:.4f}'.format(phi)+'_lamb{:.4f}'.format(lam)+'_sim'+str(sim)+'_fields.dat'
sim_suffix = lambda phi, lam, tem: '_for_temp{:.4f}'.format(temp)+'_phi0{:.4f}'.format(phi0)+'_lamb{:.4f}'.format(lamb)
bubbles_file = lambda phi, lam, tem, minsim, maxsim: pickle_file+'bubbles_from_sim'+str(minsim)+'_up_to'+str(maxsim-1)+sim_suffix(phi, lam, tem)+'.npy'
posDataFile = lambda phi, lam, tem, minsim, maxsim, mltpl, sigma: pickle_file+'positiveTargets_sims'+str(minsim)+'_to'+str(maxsim-1)+'_multiplier'+str(mltpl)+'_filter{:.4f}'.format(sigma)+sim_suffix(phi, lam, tem)+'.npy'
negDataFile = lambda phi, lam, tem, minsim, maxsim, mltpl, sigma: pickle_file+'negativeTargets_sims'+str(minsim)+'_to'+str(maxsim-1)+'_multiplier'+str(mltpl)+'_filter{:.4f}'.format(sigma)+sim_suffix(phi, lam, tem)+'.npy'
bubbles_file = lambda phi, lam, tem, minsim, maxsim: pickle_file+'bubbles_from_sim'+str(minsim)+'_up_to'+str(maxsim-1)+sim_suffix(phi, lam, tem)+'.npy'
bubble_at_rest = lambda phi, lam, tem, sim: pickle_file+'rest_bubble_sim'+str(sim)+sim_suffix(phi, lam, tem)+'.npy'
average_bubble = lambda phi, lam, tem, minsim, maxsim: pickle_file+'average_bubble_sims'+str(minsim)+'_to'+str(maxsim)+sim_suffix(phi, lam, tem)+'.npy'
full_bubble_data = lambda phi, lam, tem, minsim, maxsim: pickle_file+'rest_bubble_and_momentum_sims'+str(minsim)+'_to'+str(maxsim)+sim_suffix(phi, lam, tem)+'.npy'
average_of_N_bubbles = lambda Nbubbles, phi, lam, tem: pickle_file+'average_of_'+str(Nbubbles)+'_bubbles'+sim_suffix(phi, lam, tem)+'.npy'
VISUALISING SIMULATION DATA¶
In [5]:
def add_mask(field_slice, threshold):
return field_slice * [0 if np.abs(i) > threshold else 1 for i in field_slice]
def find_slice_peaks(field_slice, peak_threshold):
""" Finds x coordinate of peaks in masked field with mask applied at threshold. """
peak_coord = signal.find_peaks(field_slice, height = peak_threshold)[0].tolist()
if field_slice[-1] > peak_threshold and field_slice[0] > peak_threshold and field_slice[-1] != field_slice[0]:
if field_slice[0] > field_slice[-1] and field_slice[0] > field_slice[1]:
peak_coord.append(0)
elif field_slice[0] < field_slice[-1] and field_slice[-1] > field_slice[-2]:
peak_coord.append(len(field_slice)-1) # this minds potential boundary discontinuities
peak_heights = [field_slice[coord] for coord in peak_coord]
return peak_coord, field_slice.tolist().index(np.max(peak_heights))
def truncateNum(num, decimal_places):
StrNum = str(num)
p = StrNum.find(".") + 1 + decimal_places
return float(StrNum[0:p])
def time_at_fraction(bubble, frac, limit):
T, N = len(bubble), len(bubble[0])
right_phi_x = [np.sum([1 for x in slice if mask > x >= limit]) for slice in bubble]
time_list = [t if (right_phi_x[t] <= N*frac) else 0 for t in range(T)]
return next((t for t in time_list[::-1] if t != 0), 0)
def time_at_size(bubble, size, limit):
T, N = len(bubble), len(bubble[0])
right_phi_x = [np.sum([1 for x in slice if mask > x >= limit]) for slice in bubble]
time_list = [t if (right_phi_x[t] <= size) else 0 for t in range(T)]
return next((t for t in time_list[::-1] if t != 0), T-1)
def center_bubble(bubble):
limit = np.floor(phi_upper_bound)
tdecap = time_at_fraction(bubble[0], 0.85, limit)
bubble = [bubble[col][:tdecap] for col in range(len(bubble))]
bubble0 = bubble[0]
#truncate bubble where it expands relativistically
T, N = len(bubble0), len(bubble0[0])
#rotate bubble such that the relativistic bit is perfectly centered
fld = gaussian_filter1d(bubble0[-1], sigma=1, mode='wrap')
vals = [x if (fld[x] > limit and fld[x+1] > limit) else 0 for x in range(N-1)]
bubbles = [list(g) for k, g in groupby(vals, lambda x: x != 0) if k]
if len(bubbles) > 1:
maxi, imax, keep = 0, bubbles[0][0], 0
for i in range(len(bubbles)-1):
maxi = np.abs(bubbles[i+1][0]-bubbles[i][-1])
if maxi > imax:
imax = maxi
keep = i
first_zero = int(np.mean([bubbles[keep][-1], bubbles[keep+1][0]]))
bubble = np.asarray([[np.roll(slice, -first_zero) for slice in bubble[col]] for col in range(len(bubble))])
else:
first_zero = int(bubbles[0][0]/2.)
bubble = np.asarray([[np.roll(slice, -first_zero) for slice in bubble[col]] for col in range(len(bubble))])
# bubble = np.asarray([bubble[col] for col in range(len(bubble))])
# check which sides the COM travels
bubble0 = bubble[0]
tcheck = time_at_fraction(bubble0, 0.015, limit)
fld = [gaussian_filter1d(slice, sigma=1, mode='wrap') for slice in bubble0[tcheck-100:tcheck:2]]
tv = [np.nanmean([x for x in range(N) if fld[t][x] > limit]) for t in range(len(fld))]
decr = np.sum([1 for i, j in zip(tv, tv[1::]) if i > j])
incr = np.sum([1 for i, j in zip(tv, tv[1::]) if i < j])
return bubble, decr - incr #np.nanmean(tv) - N//2
def multiply_bubble(bubble, dir, fold):
bubble0 = bubble[0]
T, N = len(bubble0), len(bubble0[0])
bubble = [np.tile(bubble[col], fold) for col in range(len(bubble))] # multiplies bubbles so tail is kept without pbc
TT, NN = len(bubble[0]), len(bubble[0][0])
for t in range(TT):
a, b = int((TT-t)*light_cone) + N, int((TT-t)*light_cone/2) - N//4
x1, x2 = np.arange(a, NN), np.arange(b)
if dir < 0:
x1, x2 = x1 - a, x2 - (b-NN)
for x in np.append(x1, x2):
if 0 <= x < NN:
bubble[0][t][x] = phi_initial
return np.asarray(bubble)
In [6]:
list_temp = [0., 0.1, 0.2, 0.25, 0.3]
temp = list_temp[0]
all_data, sims_to_keep = np.load(bubbles_file(phi0, lamb, temp, lSim, nSims))
In [7]:
if False:#True:
for sim in range(len(all_data))[::-10]:
print(sim)
col = 0
simulation = all_data[sim]
centered, dir = center_bubble(simulation)
mult = multiply_bubble(centered, dir, 2)
fig, ax = plt.subplots(1, 3, figsize = (12, 3))
im0 = ax[0].imshow(simulation[col], aspect='auto', interpolation='none', origin='lower')
im1 = ax[1].imshow(centered[col], aspect='auto', interpolation='none', origin='lower')
im2 = ax[2].imshow(mult[col], aspect='auto', interpolation='none', origin='lower')
clb = plt.colorbar(im0, ax = ax[0]); clb = plt.colorbar(im1, ax = ax[1]); clb = plt.colorbar(im2, ax = ax[2]); plt.show()
In [8]:
rapidity = lambda vel: np.arctanh(vel)
gamma = lambda vel: (1.-vel**2.)**-0.5
tanh_pos = lambda x, r0L, r0R, dr, vL, vR : ( np.tanh( (x - r0L)/(dr/gamma(vL)) ) + np.tanh( (r0R - x)/(dr/gamma(vR)) ) ) * np.pi/2 + np.pi
tanh_neg = lambda x, r0L, r0R, dr, vL, vR : - ( np.tanh( (x - r0L)/(dr/gamma(vL)) ) + np.tanh( (r0R - x)/(dr/gamma(vR)) ) ) * np.pi/2 + np.pi
hyperbola1 = lambda t, a, b, c: np.sqrt(c+(t-b)**2.)+a
hyperbola2 = lambda t, d, e, f: -np.sqrt(f+(t-e)**2.)+d
def tanh_fit(bubble_slice, axis, prior):
# plt.plot(axis, bubble_slice, 'ro', axis, [tanh(r, *best_tanh) for r in axis], 'g'); plt.show()
bounds = ((axis[0], 0, 0, 0, 0), (0, axis[-1], axis[-1], 1, 1))
if prior is not None:
if np.mean(bubble_slice) > phi_initial:
return sco.curve_fit(tanh_pos, axis, bubble_slice, p0=prior, bounds=bounds)[0]
else:
return sco.curve_fit(tanh_neg, axis, bubble_slice, p0=prior, bounds=bounds)[0]
else:
if np.mean(bubble_slice) > phi_initial:
return sco.curve_fit(tanh_pos, axis, bubble_slice, bounds=bounds)[0]
else:
return sco.curve_fit(tanh_neg, axis, bubble_slice, bounds=bounds)[0]
def hypfit_test(tt, rr, qq): # no parameters in common between left and right wall for hyperbolic curve fit
try:
if rr[0] <= rr[-1]:
fit1, pcov1 = sco.curve_fit(hyperbola1, tt, rr, p0 = (min(rr), tt[rr.tolist().index(min(rr))], 1e4))
fit2, pcov2 = sco.curve_fit(hyperbola2, tt, qq, p0 = (max(qq), tt[qq.tolist().index(max(qq))], 1e4))
else:
fit1, pcov1 = sco.curve_fit(hyperbola2, tt, rr, p0 = (max(rr), tt[rr.tolist().index(max(rr))], 1e4))
fit2, pcov2 = sco.curve_fit(hyperbola1, tt, qq, p0 = (min(qq), tt[qq.tolist().index(min(qq))], 1e4))
return True
except (RuntimeError, TypeError):
return False
def hypfit(tt, rr, qq): # no parameters in common between left and right wall for hyperbolic curve fit
if rr[0] <= rr[-1]:
fit1, pcov1 = sco.curve_fit(hyperbola1, tt, rr, p0 = (min(rr), tt[rr.tolist().index(min(rr))], 1e4))
fit2, pcov2 = sco.curve_fit(hyperbola2, tt, qq, p0 = (max(qq), tt[qq.tolist().index(max(qq))], 1e4))
traj1, traj2 = hyperbola1(tt, *fit1), hyperbola2(tt, *fit2)
else:
fit1, pcov1 = sco.curve_fit(hyperbola2, tt, rr, p0 = (max(rr), tt[rr.tolist().index(max(rr))], 1e4))
fit2, pcov2 = sco.curve_fit(hyperbola1, tt, qq, p0 = (min(qq), tt[qq.tolist().index(min(qq))], 1e4))
traj1, traj2 = hyperbola2(tt, *fit1), hyperbola1(tt, *fit2)
# print('hyperbolic trajectories fit: ', fit1, fit2)
# plt.plot(rr, tt, 'g-', qq, tt, 'r-', traj1, tt, 'y:', traj2, tt, 'b:') # plot the equation using the fitted parameters
return tt, traj1, traj2
def get_velocities(tt, rrwallfit, llwallfit):
#uu = wall travelling along with the com; vv = wall travelling against the com; aa = centre of mass velocity; bb = instantaneous velocity of wall
uu = intp.splev(tt, intp.splrep(tt, rrwallfit), der=1)
vv = intp.splev(tt, intp.splrep(tt, llwallfit), der=1)
for i in range(len(uu)):
if np.abs(uu[i]) > 1:
uu[i] = np.sign(uu[i])*(1-np.abs(np.abs(uu[i])-1))
if np.abs(vv[i]) > 1:
vv[i] = np.sign(vv[i])*(1-np.abs(np.abs(vv[i])-1))
# print(uu)
# print(vv)
for i in range(len(uu)):
if str(uu[i]) == 'nan' :
uu[i] = - np.sign(vv[-1])*0.999
if str(vv[i]) == 'nan' :
vv[i] = - np.sign(uu[-1])*0.999
aa = ( 1 + uu*vv - np.sqrt( (-1 + uu**2)*(-1 + vv**2))) / ( uu + vv)
bb = (-1 + uu*vv + np.sqrt( (-1 + uu**2)*(-1 + vv**2))) / (-uu + vv)
return uu, vv, aa, bb#, ch1, ch2
def velocity(bubble, bool1, bool2):
T, N = len(bubble), len(bubble[0])
limit = phi_upper_bound #phi_upper_bound - 0.95*np.abs(phi_upper_bound-np.floor(phi_upper_bound))
err = 0
# find where fit can begin
window = int(N//10)
if window > 100: window = 100
tf = T-1
endSlice = [i for i in range(N) if bubble[tf][i] > limit]
if len(endSlice) == 0:
return 'next'
while endSlice[0] - window < 0 or endSlice[-1] + window >= N:
tf -= 2
endSlice = [i for i in range(N) if bubble[tf][i] > limit]
if len(endSlice) == 0:
return 'next'
data_list = []
prior = None
for t in range(tf, 0, -1):
if len(data_list) == 0:
peaks = endSlice
target = int(np.round(np.nanmean(peaks)))
coord_list = np.arange(peaks[0] - window, peaks[-1] + window)
else:
xrange = np.arange(data_list[-1][0]-window, data_list[-1][0]+window+1)
if all(0 <= i < N for i in xrange) and any(np.cos(bubble[t][i]) > np.cos(limit) for i in xrange):
peaks = [i for i in xrange if np.cos(bubble[t][i]) > np.cos(limit)]
else:
break
target = int(np.round(np.nanmean(peaks)))
coord_list = np.arange(target - int(np.round(np.abs(data_list[-1][1]))) - window, target + int(np.round(np.abs(data_list[-1][2]))) + window)
coords = np.asarray([(bubble[t][i%N], int(i-target)) for i in coord_list])
if err < 5:
try:
r0L, r0R, dr, vL, vR = tanh_fit(coords[:,0], coords[:,1], prior)
prior = None#r0L, r0R, dr, vL, vR
data_list.append([target, r0L, r0R, int(t)])
except (RuntimeError, ValueError, TypeError):
prior = None
data_list.append([target, r0L + np.sign(data_list[-1][1]-data_list[0][1])*light_cone*(T/N), r0R + np.sign(data_list[-1][2]-data_list[0][2])*light_cone*(T/N), int(t)])
err += 1
else:
break
data_list = np.asarray(data_list[::-1])
targets, r0Ls, r0Rs, time_list = data_list[:,0], data_list[:,1], data_list[:,2], data_list[:,-1]
# get direction of bubble
radius_diff = np.mean([np.abs(data_list[i,1]) - np.abs(data_list[i,2]) for i in range(len(data_list)//4)])
if radius_diff > 0: # if average difference in radius is positive then left radius is on average higher than right radius i.e. bubble travels to the right
rr, ll = targets + r0Rs, targets + r0Ls
else:
rr, ll = targets + r0Ls, targets + r0Rs
if bool1:
fig, ax0 = plt.subplots(1, 1, figsize = (5, 4))
ax0.plot(rr, time_list, 'b', ll, time_list, 'y', linewidth='3')
ax0.set_xlabel('x'); ax0.set_ylabel('t'); plt.show()
# get velocities from derivative of best fit to wall trajectory
trunc = 0
# save copies
time_list_copy, rr_copy, ll_copy = time_list, rr, ll
# as lons as it needs..
while(True):
# try to fit walls
if hypfit_test(time_list_copy, rr_copy, ll_copy):
time_list, rrwallfit, llwallfit = hypfit(time_list_copy, rr_copy, ll_copy)
uu,vv,aa,bb = get_velocities(time_list, rrwallfit, llwallfit)
break
# if fit fails, shorten walls
else:
trunc += 25
if len(time_list_copy) > 100:
time_list_copy, rr_copy, ll_copy = time_list_copy[trunc::], rr_copy[trunc::], ll_copy[trunc::]
else:
return 'next'
continue
if bool2:
fig, [ax0, ax1] = plt.subplots(1, 2, figsize = (15, 4))
ax0.plot(rr[-len(time_list):], time_list, 'b', ll[-len(time_list):], time_list, 'y', rrwallfit, time_list, 'r:', llwallfit, time_list, 'r:', linewidth='3')
ax0.set_xlabel('x'); ax0.set_ylabel('t')
ax1.plot(time_list, uu, 'b:', time_list, vv, 'y:')
ax1.plot(time_list, aa, 'r', label=r'v COM')
ax1.plot(time_list, bb, 'g', label=r'v walls')
ax1.axhline(0, color='darkgray', ls=':')
ax1.set_xlabel('t'); ax1.set_ylabel('v(t)'); ax1.legend(); plt.show()
list = np.abs(aa - bb)
if any(str(i) != 'nan' for i in list):
return -aa[list.tolist().index(np.nanmin(list))]
else:
return 'next'
In [9]:
fold = lambda beta: 3 if np.abs(beta) > 0.8 else 2
add_velocities = lambda v1, v2: (v1 + v2) / (1. + v1*v2)
def cut_out_bubble(bubble):
bubble0 = bubble[0]
limit = np.floor(phi_upper_bound)
T, N = len(bubble0), len(bubble0[0])
ext = int(N//8)
if ext < 100: ext = 100
xmin = next((i for i in range(N) if bubble0[-1][i] > limit), 0) - ext
xmax = next((i for i in range(N)[::-1] if bubble0[-1][i] > limit), N-1) + ext
if xmin < 0: xmin = 0
if xmax >= N: xmax = N-1
t, tmin = 0, 0
while not any(bubble0[t][i] > limit for i in range(xmin, xmax+1)):
tmin = t
t += 1
return [[[bubble[col][t][x] for x in range(xmin, xmax+1)] for t in range(tmin, T)] for col in range(len(bubble))]
def center_bubble_again(bubble, kind):
bubble0 = bubble[0]
T, N = len(bubble0), len(bubble0[0])
ext = int(N//8)
if ext < 100: ext = 100
if kind < 0:
x1 = 0
peaks1 = [phi_initial, phi_initial]
while np.cos(np.nanmean(peaks1)) < np.cos(3.18) and x1 < N-1:
# while len([i for i,j in zip(peaks1, peaks1[1:]) if (i > right_phi_at_V_max.x and j > right_phi_at_V_max.x)]) > 10 and x1 < N-1:
slice = [bubble0[i][x1] for i in range(T)]
peaks1 = [i for i in slice if i < mask]
peaks2 = [i for i in range(T) if slice[i] < mask]
x1 += 1
if len(peaks2) > 0:
tmax = peaks2[-1]
else:
return 'next'
# print('tmax, x1', tmax, x1)
xmin = x1 - ext
if xmin < 0: xmin = 0
xmax = next((i for i in range(N-1, xmin, -1) if (mask > bubble0[tmax][i] > phi_initial)), N-1)
else:
x1 = N-1
peaks1 = [phi_initial, phi_initial]
while np.cos(np.nanmean(peaks1)) < np.cos(3.18) and x1 > 0:
# while len([i for i,j in zip(peaks1, peaks1[1:]) if (i > right_phi_at_V_max.x and j > right_phi_at_V_max.x)]) > 10 and x1 > 0:
slice = [bubble0[i][x1] for i in range(T)]
peaks1 = [i for i in slice if i < mask]
peaks2 = [i for i in range(T) if slice[i] < mask]
x1 -= 1
if len(peaks2) > 0:
tmax = peaks2[-1]
else:
return 'next'
# print('tmax, x1', tmax, x1)
xmax = x1 + ext
if xmax > N-1: xmax = N-1
xmin = next((i for i in range(xmax) if (mask > bubble0[tmax][i] > phi_initial)), 0)
t, tmin = 0, 0
while not any([(mask > bubble0[t][i] > phi_initial) for i in range(xmin, xmax+1)]):
tmin = t
t += 1
return [[[bubble[col][t][x] if bubble[col][t][x] != mask else normal[col] for x in range(xmin, xmax+1)] for t in range(tmin, tmax+1)] for col in range(len(bubble))]
def interpolate_bubble(bubble, res):
NT, N = len(bubble), len(bubble[0])
bubble = np.asarray(bubble)
t, x = np.arange(NT), np.arange(N)
tt, xx = np.meshgrid(t, x)
f = intp.interp2d(t, x, bubble[tt, xx], kind = 'quintic')
T, X = np.arange(0, NT, 1/res), np.arange(0, N, 1/res)
return f(T, X).T
def boost_bubble(bubble, vCOM, res):
print('Boost velocity factor: ', vCOM)
boost = lambda vCOM, ga, x, t: ga * (x - vCOM*t)
bubble0 = bubble[0]
T, N = len(bubble0), len(bubble0[0])
limit = phi_upper_bound
ga = gamma(vCOM)
t0 = int(time_at_size(bubble0, 50, limit)) # arbitrary choice of bubble centre
x0 = int(np.mean([i for i in range(N) if mask > bubble0[t0][i] > limit]))
print('boost centre t0, x0: ', t0, x0)
T0, X0 = boost(vCOM, ga, t0, x0), boost(vCOM, ga, x0, t0)
deltaT, deltaX = T0-t0, X0-x0
bubble = [interpolate_bubble(bubble[col], res) for col in range(len(bubble))]
new_bubble = [[[mask for x in range(N)] for t in range(T)] for col in range(len(bubble))]
for t in range(T):
for x in range(N):
tt = int((boost(vCOM, ga, t, x) + np.sign(T/2.-t0) * deltaT)*res)
xx = int((boost(vCOM, ga, x, t) + np.sign(N/2.-x0) * deltaX)*res)
if (T*res > tt >= 0 and N*res > xx >= 0):
for col in range(len(bubble)):
new_bubble[col][t][x] = bubble[col][tt][xx]
return new_bubble, x0 - int(N/2.)
In [10]:
def the_whole_shebang(resolution, saveF):
for sim in range(len(all_data))[::]:
exists = os.path.exists(bubble_at_rest(phi0, lamb, temp, sims_to_keep[sim]))
if sims_to_keep[sim] >= 0 and not exists:
#if True:
print('On to simulation ', sim)
bubble = all_data[sim]
bubble, dir = center_bubble(bubble)
beta = velocity(bubble[0], bool1=False, bool2=True)
if beta == 'next':
print('sim '+str(sim)+' skipped')
continue
elif np.abs(beta) > 0.9:
beta = np.sign(beta)*0.9
bubble = multiply_bubble(bubble, dir, fold(beta))
fig = plt.figure(figsize=(5,5))
plt.imshow(bubble[0], aspect='auto', interpolation='none', origin='lower')
plt.title('original view'); plt.show()
vtotal, bool1, bool2 = 0, True, False
while np.abs(beta) > 0.1:
bool2 = True
bubble, dir = boost_bubble(bubble, beta, resolution)
fig = plt.figure(figsize=(5,5))
plt.imshow(bubble[0], aspect='auto', interpolation='none', origin='lower')
plt.title('boosted view'); plt.show()
bubble = center_bubble_again(bubble, dir)
if bubble == 'next' or any(np.shape(bubble[0])) == 0:
print('sim '+str(sim)+' interrupted')
bool1 = False
break
else:
fig = plt.figure(figsize=(5,5))
plt.imshow(bubble[0], aspect='auto', interpolation='none', origin='lower')
plt.title('re-centered view'); plt.show()
prebeta = beta
vtotal = add_velocities(vtotal, beta)
print('step done ', beta)
min_bubble = cut_out_bubble(bubble)
fig = plt.figure(figsize=(5,5))
plt.imshow(min_bubble[0], aspect='auto', interpolation='none', origin='lower')
plt.title('final view'); plt.show()
beta = velocity(min_bubble[0], bool1=False, bool2=True)
if beta == 'next':
print('sim '+str(sim)+' interrupted')
if np.abs(prebeta) >= 0.8:
bool1 = False
break
elif np.abs(beta) > 0.9:
beta = np.sign(beta)*0.9
if not bool2:
if np.abs(beta) > 0:
min_bubble = cut_out_bubble(bubble)
vtotal = beta
bool2 = True
if bool1 and bool2:
print('sim', sim, 'total vel ', vtotal)
if saveF:
np.save(bubble_at_rest(phi0, lamb, temp, sims_to_keep[sim]), [min_bubble, vtotal])
print('Bubble saved!')
return
In [11]:
#the_whole_shebang(2, saveF = True)
the_whole_shebang(2, saveF = True)
Average bubble¶
In [12]:
def cut_out_coordinates(bubble0):
T, N = len(bubble0), len(bubble0[0])
limit = phi_upper_bound
size = 50 # do not change
tmin = T-1
for t in range(T):
slice = gaussian_filter1d(bubble0[t], sigma=5, mode='wrap')
right_phi = [x for x in range(N-1) if (slice[x] >= limit and slice[x+1] >= limit)]
if len(right_phi) > size:
tmin = t
break
if tmin == T-1:
return None
xcen = int(np.round(np.nanmean(right_phi)))
xcen_list = [np.round(np.nanmean([x for x in range(N) if bubble0[t][x] >= limit])) for t in range(tmin-50,tmin,2)]
if stat.stdev(xcen_list) > 90:
return None
# fig = plt.figure(figsize=(4,4))
# plt.plot(xcen, tmin, 'ro')
# plt.imshow(bubble, aspect='auto', interpolation='none', origin='lower')
# plt.show()
ex1, ex2, ex3, ex4 = np.arange(tmin)[::-1], np.arange(tmin, T), np.arange(xcen)[::-1], np.arange(xcen, N)
if len(ex1) > 1500: ex1 = ex1[:1500]
if len(ex2) > 500: ex2 = ex2[:500]
if len(ex3) > 500: ex3 = ex3[:500]
if len(ex4) > 500: ex4 = ex4[:500]
return ex1, ex2, ex3, ex4
def stack_bubbles(bubble_data, bool1, bool2):
extent_list, bad_bubbles = [], []
for sim in range(len(bubble_data)):
coords = cut_out_coordinates(bubble_data[sim][0])
if coords is not None:
extent_list.append(coords)
else:
bad_bubbles.append(sim)
if bool1:
fig, ax = plt.subplots(1, 1, figsize=(6, 5))
im = ax.imshow(bubble_data[sim][0], aspect='auto', interpolation='none', origin='lower')
clb = plt.colorbar(im, ax = ax); plt.show()
extent_list = np.asarray(extent_list)
print('bad bubbles', bad_bubbles, ' pictured above if bool1 = True')
bubble_data = [bubble_data[i] for i in range(len(bubble_data)) if i not in bad_bubbles]
print('Number of good bubbles remaining: ', len(bubble_data), ' pictured below if bool2 = True')
#plot good bubbles remaining
if bool2:
for sim in range(len(bubble_data)):
fig, ax = plt.subplots(1, 1, figsize=(6, 5))
# print(extent_list[sim][0][-1], extent_list[sim][1][-1], extent_list[sim][2][-1], extent_list[sim][3][-1])
im = ax.imshow(np.asarray(bubble_data[sim][0])[extent_list[sim][0][-1]:extent_list[sim][1][-1], extent_list[sim][2][-1]:extent_list[sim][3][-1]], aspect='auto', interpolation='none', origin='lower')
clb = plt.colorbar(im, ax = ax); plt.show()
# for sim in range(len(bubble_data)):
# bubble_data[sim].append([0.5*np.asarray(i)**2 for i in bubble_data[sim][1]])
nCols = len(bubble_data[0])
uper_right_bubble_list = [[[[bubble_data[sim][col][t][x] for t in extent_list[:,0][sim]] for x in extent_list[:,2][sim]] for sim in range(len(bubble_data))] for col in range(nCols)]
upper_left_bubble_list = [[[[bubble_data[sim][col][t][x] for t in extent_list[:,0][sim]] for x in extent_list[:,3][sim]] for sim in range(len(bubble_data))] for col in range(nCols)]
lower_right_bubble_list = [[[[bubble_data[sim][col][t][x] for t in extent_list[:,1][sim]] for x in extent_list[:,2][sim]] for sim in range(len(bubble_data))] for col in range(nCols)]
lower_left_bubble_list = [[[[bubble_data[sim][col][t][x] for t in extent_list[:,1][sim]] for x in extent_list[:,3][sim]] for sim in range(len(bubble_data))] for col in range(nCols)]
return [uper_right_bubble_list, upper_left_bubble_list, lower_right_bubble_list, lower_left_bubble_list], bad_bubbles
def average_bubble_func(stacked_bubble_data):
uper_right_bubble_list, upper_left_bubble_list, lower_right_bubble_list, lower_left_bubble_list = stacked_bubble_data
nCols = len(uper_right_bubble_list)
print(nCols)
av_mat, av_err_mat = [], []
for col in range(nCols):
average_matrix, average_error_matrix = [[],[],[],[]], [[],[],[],[]]
coord_lists = [uper_right_bubble_list[col], upper_left_bubble_list[col], lower_right_bubble_list[col], lower_left_bubble_list[col]]
for i in range(4): #for each quadrant
max_lines, max_cols = 0, 0
for simulation in coord_lists[i]:
if len(simulation) > max_lines:
max_lines = len(simulation)
if len(simulation[0]) > max_cols:
max_cols = len(simulation[0])
for line in range(max_lines):
average_matrix[i].append([normal[col] for nn in range(max_cols)])
average_error_matrix[i].append([0. for nn in range(max_cols)])
for num_line in range(max_lines):
for num_col in range(max_cols):
meas = []
for simulation in coord_lists[i]:
if len(simulation) > num_line and len(simulation[0]) > num_col:
meas.append(simulation[num_line][num_col])
average_matrix[i][num_line][num_col] = np.mean(meas)
if len(meas) >= 2:
average_error_matrix[i][num_line][num_col] = stat.stdev(meas)/len(meas)
av_mat.append(average_matrix)
av_err_mat.append(average_error_matrix)
whole_bubbles = []
for av_bub in [av_mat, av_err_mat]:
whole_bubbles.append([])
for col in range(nCols):
top = np.concatenate((np.flip(np.asarray(np.flip(av_bub[col][0],1)).transpose(),1), np.flip(np.asarray(av_bub[col][1]).transpose(),0)), axis=1)
bottom = np.concatenate((np.flip(np.asarray(av_bub[col][2]).transpose(),1), np.flip(np.asarray(np.flip(av_bub[col][3],0)).transpose(),1)), axis=1)
whole_bubbles[-1].append(np.concatenate((top, bottom), axis=0))
return np.asarray(whole_bubbles)
def truncated_bubbles_stacked(stacked_bubble_data):
uper_right_bubble_list, upper_left_bubble_list, lower_right_bubble_list, lower_left_bubble_list = stacked_bubble_data
nCols = len(uper_right_bubble_list)
whole_bubbles = []
for sim in range(len(uper_right_bubble_list[0])):
whole_bubbles.append([])
for col in range(nCols):
bub0 = uper_right_bubble_list[col][sim]
bub1 = upper_left_bubble_list[col][sim]
bub2 = lower_right_bubble_list[col][sim]
bub3 = lower_left_bubble_list[col][sim]
top = np.concatenate((np.flip(np.asarray(np.flip(bub0,1)).transpose(),1), np.flip(np.asarray(bub1).transpose(),0)), axis=1)
bottom = np.concatenate((np.flip(np.asarray(bub2).transpose(),1), np.flip(np.asarray(np.flip(bub3,0)).transpose(),1)), axis=1)
whole_bubbles[-1].append(np.concatenate((top, bottom), axis=0))
return whole_bubbles
In [13]:
bubble_list = np.asarray([np.load(bubble_at_rest(phi0, lamb, temp, sim)) for sim in range(lSim, nSims) if os.path.exists(bubble_at_rest(phi0, lamb, temp, sim))])
print('Working with ', len(bubble_list), ' bubbles at rest.')
bubble_list_save, bubble_velocities = bubble_list[:,0], bubble_list[:,1]
del bubble_list
In [14]:
#either compute average
stacked_bubble_parts, bad_bubbles = stack_bubbles(bubble_list_save, bool1=True, bool2=True)
# save data
ab, error_ab = average_bubble_func(stacked_bubble_parts)
newBubbleData = truncated_bubbles_stacked(stacked_bubble_parts)
np.save(average_bubble(phi0, lamb, temp, lSim, nSims), [ab, error_ab])
np.save(full_bubble_data(phi0, lamb, temp, lSim, nSims), newBubbleData)
# or load it
#ab, error_ab = np.load(average_bubble(phi0, lamb, temp, lSim, nSims))
#newBubbleData = np.load(full_bubble_data(phi0, lamb, temp, lSim, nSims)).tolist()
#separate columns
bubble = ab#, mom_bubble, ge_bubble, pe_bubble, ke_bubble = ab
error_bubble = error_ab#, error_mom_bubble, error_ge_bubble, error_pe_bubble, error_ke_bubble = error_ab
In [15]:
#plot average bubble, mom_bubble, ge_bubble, pe_bubble, ke_bubble
fig, ax = plt.subplots(1, len(ab)+1, figsize=(12, 5))
for col in range(len(ab)):
im = ax[col].imshow(ab[col], aspect='auto', interpolation='none', origin='lower')
clb = plt.colorbar(im, ax = ax[col])
plt.savefig(plots_file+'average_bubble'+sim_suffix(phi0, lamb, temp)+'.png')
#plot distribution of boosts
bubble_velocities = [bubble_velocities[i] for i in range(len(bubble_velocities)) if i not in bad_bubbles]
fig, ax = plt.subplots(1, 3, figsize=(15,5))
ax[0].hist(bubble_velocities, bins=len(bubble_velocities))
ax[1].hist([gamma(i) for i in bubble_velocities], bins=len(bubble_velocities))
ax[2].hist([rapidity(i) for i in bubble_velocities], bins=len(bubble_velocities))
[ax[i].set_ylabel('Count') for i in range(len(ax))]
ax[0].set_xlabel(r'$v$'); ax[1].set_xlabel(r'$\gamma$'); ax[2].set_xlabel(r'$w$');
plt.savefig(plots_file+'boost_distrib'+sim_suffix(phi0, lamb, temp)+'.png');
In [16]:
#average x bubbles at a time; save
Nbubs = len(newBubbleData)
step = 2
for ii in range(0, Nbubs, step):
boolBreak = False
qq = len(newBubbleData)
print(len(newBubbleData))
stacked_bubble_parts, _ = stack_bubbles(newBubbleData, False, False)#np.asarray(newBubbleData)[:,:-1].tolist())
print(len(newBubbleData))
ab, error_ab = average_bubble_func(stacked_bubble_parts)
np.save(average_of_N_bubbles(qq, phi0, lamb, temp), [ab, error_ab, qq])
for jj in range(step):
if len(newBubbleData) > 2:
del newBubbleData[0]
else:
boolBreak = True
break
if boolBreak:
break
In [17]:
image1 = np.asarray(ab[0])[800:1100, 300:700]
fig, ax = plt.subplots(1, 1, figsize = (5, 3))
im = ax.imshow(image1, aspect='auto', interpolation='none', origin='lower')
plt.show()
beta = velocity(image1, bool1 = False, bool2 = True)
In [18]:
fig, ax = plt.subplots(1, 2, figsize = (12,4), gridspec_kw={'width_ratios': [4, 5]})
liml = 160
limm = 200
ext = [0, image1.shape[1], 0, image1.shape[0]]
#ext = [0, image1.shape[1], 0, image1.shape[0]]
im = ax[0].imshow(image1, aspect='auto', interpolation='none', origin='lower', cmap='viridis', extent=ext)
ax[0].axhline(liml, color='r', linestyle='--')
ax[0].axhline(limm, color='r', linestyle='--')
clb = plt.colorbar(im, ax = ax[0], pad=0.01); clb.set_label(r'$\bar{\phi}$')
multp = 3
clb.set_ticks(np.arange(10)*np.pi/multp, update_ticks=True)
clb.set_ticklabels([r'$\pi$' if (iii/multp)==1. else str(int(iii/multp))+r'$\pi$' if (iii/multp).is_integer()
else r'$\frac{%s\pi}{%s}$'%(iii, multp) for iii in np.arange(10)], update_ticks=True)
ax[0].set(xlabel = r'$\phi_0^{-1} \sqrt{V_0} \; r$', ylabel = r'$\phi_0^{-1} \sqrt{V_0} \; t$')
saves = []
for temp in image1[liml:limm]:
ax[1].plot(np.arange(len(temp)), temp, color='b', alpha=0.05, linewidth=1)
peaks,_ = scs.find_peaks(temp)
fwhm, height, left_ips, right_ips = scs.peak_widths(temp, peaks, rel_height=0.5)
# fwhm, left_ips, right_ips, peaks = fwhm[height==max(height)], left_ips[height==max(height)], right_ips[height==max(height)], peaks[height==max(height)]
# saves.append([temp, fwhm/2., max(height), left_ips, right_ips, peaks])
height, left_ips, right_ips, peaks = height[fwhm==max(fwhm)], left_ips[fwhm==max(fwhm)], right_ips[fwhm==max(fwhm)], peaks[fwhm==max(fwhm)]
ax[1].plot(peaks, temp[peaks], marker="o", color='green', ms=4)
saves.append([temp, max(fwhm)/2., height, left_ips, right_ips, peaks])
saves = np.asarray(saves)
temps, fwhms, heights, lefts, rights, peaks = saves[:,0], saves[:,1], saves[:,2], saves[:,3], saves[:,4], saves[:,5]
ax[1].plot(np.arange(len(np.mean(temps))), np.mean(temps), color='r', alpha=1)
pk, _ = scs.find_peaks(np.mean(temps))
w, h, l, r = scs.peak_widths(np.mean(temps), pk, rel_height=0.5)
h, l, r, pk = h[w==max(w)], l[w==max(w)], r[w==max(w)], pk[w==max(w)]
#l, r, pk = l[h==max(h)], r[h==max(h)], pk[h==max(h)]
ax[1].hlines(h, l, r, color='orange', linewidth=3)
ax[1].plot(pk, np.mean(temps)[pk], marker="*", color='green', ms=12)
#ax[1].hlines(np.mean(heights), np.mean(lefts), np.mean(rights), color=tableau1[0], linewidth=3)
#ax[1].plot(np.mean(peaks), np.mean(temps)[np.mean(peaks)], marker="*", color=tableau1[2], ms=12)
plt.grid(alpha=0.8, linestyle='dashed', linewidth=0.5)
ax[1].set_xlabel(r'$\phi_0^{-1} \sqrt{V_0} \; r$'); ax[1].set_ylabel(r'$\bar{\phi}$')
plt.savefig(plots_file+'average_thermal_bubble.pdf');
plt.show()
fwhms = fwhms[(np.mean(fwhms)-4*np.std(fwhms) < fwhms) & (fwhms < np.mean(fwhms)+4*np.std(fwhms))]
filter_size = np.mean(fwhms); print('filter size = ', filter_size)
smoothen = lambda slice, sigma: np.fft.ifft(np.fft.fft(slice)* np.exp(-0.5* (np.fft.fftfreq(len(slice), dx)*2*np.pi * dx * sigma)**2.)).real
MAXlist = [max(smoothen(temp, filter_size)) for temp in image1[liml:limm]]
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: